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1. INTRODUCTION 

Two main controlling factors of ice crystal growth 
are the heat and (gaseous) mass transfers, H&MT, 
as characterized by the Nusselt number (Nu) for heat 
and the Sherwood (Sh) number for mass. Nu and 
Sh express the increase of the molecular diffusions by 
the relative motions of the particles in air as charac- 
terized by the Reynolds number (Re). Traditionally 
Nu is assumed to be identical to Sh. Past labora- 
tory measurements with various ice crystal shapes 
by the Toronto Cloud Physics Group involved Sh, 
thus this study will deal with Sh. The growth of a 
crystal depends on its size and shape, its fall speed 
(through Re) and the flux of water vapour from the 
surrounding air. It is controlled by the diffusion of 
heat which carries away the energy released on the 
crystal surface by deposition. Sh is a function of 
Re and was established by Schemenauer and List 
(1978)[S&L78] for snow crystals and graupel. Their 
study, however, did not address the case of pure dif- 
fusion (i.e. Re=0). Such values based on approxima- 
tions have been available in the literature (McDon- 
ald, 1963[Mc63], Jayaweera (1971)) only for a very 
limited number of shapes of crystals. Thus, it was 
decided to take a general approach to numerically 
calculate Sh^ e= o for any ice crystal with a recti- 
linear shape. New values of ShR e= o for additional 
shapes of crystals are addressed in this study along 
with the method developed for the computations. 

The steady state diffusion is controlled by the 
Laplace equation. It was solved numerically with 
Dirichlet boundary conditions on a rectilinear 3 di- 
mensional lattice system with variable lattice separa- 
tion distances for hexagonal plates, hexagonal cylin- 
ders, stellar crystals, capped columns, and broad- 
branched crystals. 



2. ELECTROSTATIC ANALOGY 

The water vapour density field (p) around a sta- 
tionary ice crystal growing under steady state diffu- 
sion obeys the Laplace equation 



V 2 p = 



(1) 



The initial growth of an ice crystal by the net de- 
position of water vapour in steady state diffusion, is 
controlled by the rate at which vapour diffuses to 
and from the crystal's surface. The rate at which 
heat is released by the crystal is exactly the same as 
the rate at which heat is diffused from the surface. 
This balance leads to quasi constant value of crystal's 
surface temperature, and this in turn leads to a con- 
stant local vapour density po around the crystal. For 
a perfect capacitor that is geometrically similar to the 
crystal, the electrostatic potential V in vacuum also 
obeys the Laplace equation, with the surface electro- 
static potential Vq constant as is the case for perfect 
capacitors. Identifying V with p, the direct analogy 
between steady state diffusion and electrostatics is 
established. A capacitor that is geometrically similar 
to the crystal is then set in a rectangular box, a Fara- 
day cage, sufficiently afar from its edges. The capac- 
itor is then surrounded by the electrostatic potential 
obeying Laplace equation. The sides of the cage are 
set at potential zero to simulate "boundaries at infin- 
ity" . Since the water vapour tapers off to a constant 
value sufficiently far away from the crystal, imposing 
this Dirichlet boundary condition preserves the sim- 
ilarity. In this analogy, one uses Gauss' flux law and 
Fick's law of diffusion to establish that Sh,R e= o is 
given by 

„, 4nCL 

A/lRe=0 — — — — (2) 

where A, L, and C are the surface area, charac- 
teristic length, and the crystal's capacitance, respec- 



tively. The units are in cgs for ease of comparison 
with S&L78 and Jayaweera(1971). Thus the prob- 
lem of diffusion is reduced to finding capacitances of 
crystals with the Dirichlet boundary conditions men- 
tioned above. 

3. FINITE LATTICES AND JACOBI'S 
METHOD 

To compute (2), a numerical scheme for finding 
the capacitances for crystals of rectilinear shapes was 
developed. "Rectilinear shapes" are those that can be 
constructed with a finite number of straight edges. 

A Cartesian grid system is established with a finite 
number of lattice points within the Faraday cage, the 
outer most points representing the edges of the cage 
(box). Each side of the box is at least 20 times larger 
than the respective sides of the crystal. The separa- 
tion between adjacent lattice points can be different 
along different axis directions. Each of the lattice 
points are labelled in terms of positive integer in- 
dices (i, j, k); where i, j, and k run along the x-axis, 
y-axis, and z-axis respectively (Fig. 1). For sim- 
plicity, ji is introduced, which can be either x, y, or 
z. Then fi represents the unit vector along /i-axis. 
Assuming that the grid points are sufficiently close 
to one another, the partial derivatives of potential 
V can be approximated by the following set of finite 
difference equations 



dV(i,j,k) _ V((i,j, k) + /}) - V(i,j, k) 



(3) 



where 



A{f£) ee V((i,j, k) + /x) + V((i,j, k) - A) (4) 

then the finite difference approximation of 2nd order 
partial derivatives of V can be written as 



d 2 V(i,j,k) A(»)-2V(i,j,k) 



2(^)2 



(5) 



For convenience, following proportionality constants 
are introduced: 



■ s <!> 2 



Then the Laplace's equation in our discrete rectilin- 
ear grid system becomes 

n ~ 3eL + B ^ + B ^ (?) 

~ 2(Sx) 2 2a(Sx) 2 2/?(<5x) 2 1 ' 



where B(/i)= A(^) - 2V(i,j,k). From (7), one can 
solve for potential V(i,j, k) at any lattice point 
(i,j,k) in our Faraday cage. By rearranging the terms 
in (7), the following is obtained: 

_ a(3A(x)+pA(y)+aA(z) 

This is the expression used by the "Ja- 
cobi's" iterative scheme for solving Laplace equation. 
An algorithm is constructed that runs with (8). First, 
the values of V are assigned to the sides of the Fara- 
day cage and the surface of the crystal, which is 
sitting in the middle of this Faraday cage. V is then 
set to zero on all other lattice points, followed by 
the running of a computational scheme for V(i,j,k) 
iteratively for each lattice point in the box at each 
iteration. The computation is started out from the 
crystal surface, then "branches out towards the edges 
of the box". This procedure is iterated until con- 
vergence by computing the difference 5V(i,j,k) — 
\V n (i,j, k) — V n -i(i,j, k)\, where V n (i, j, k) denotes 
the value V(i,j,k) at n-th iteration. The scheme is 
halted when a desired "tolerance bound" is reached. 
This occurs when following criterion is met 



E 5V ^ 



j, k) < 



(9) 



where the sum in (9) is carried out for all grid points 
in the box. It is well known from literatures in com- 
putational physics that the number of iterations V 
required to reduce the error given by the sum (9) by 
a factor of 1(T P is 



r - 0(pN 3 



(10) 



where N is the number of grid points in a NxNxN 
cubic regular Cartesian grid system. At the end of 
a run of this algorithm, values of all the V's would 
have been computed for each lattice point up to the 
desired precision e. 

Some comments are in order at this stage. Al- 
though no regular Cartesian grid system is used for 
many of the crystal models studied here, the order 
of convergence will be the same as (10). There 
are other computational methods to solve Laplace's 
equation numerically, such as Gauss-Seidel and the 
"simultaneous over-relaxation" (SOR) methods. The 
SOR method has a convergence rate of 0(N 2 ) for 
3 dimensional Cartesian grid system. The classical 



Jacobi method, however, is one of the simplest al- 
gorithms that can be used for our study and is rel- 
atively free of the complications in implementation 
compared to the SOR method. It was thus selected 
for this study. 

4. EXAMPLE: HEXAGONAL PLATE 

The hexagonal plate was modelled using a non- 
regular rectangular grid system, with the proportion- 
ality factors a and (3 not being equal. Only a fi- 
nite number of grid points is available for "drawing" a 
hexagon to encapsulate the important features of the 
crystal, namely its vertices. Not all the vertices of 
the hexagon lie on lattice points if a regular (a = (3) 
Cartesian grid system is used. Hence a regular Carte- 
sian coordinate system is not used for an optimal 
design of the crystal on our finite grid. a=3 allows 
for construction of a hexagon in which all its vertices 
lie on grid points. (3 can be chosen to be any pos- 
itive value since the plate is assumed to be thinner 
than the lattice separation along the z-axis. Thus 
(3=1 was selected, meaning that our plate lies on a 
single xy-planar layer of lattice points, in the mid- 
dle of the box. This "optimal" representation of the 

N x 

hexagonal plate results in the relationship: Nd = 

where Nd and N x represent the number of points 
used to represent a diagonal segment, and horizon- 
tal segment respectively. In a similar manner, "op- 
timal" representation was achieved for other crystal 
shapes as investigated experimentally by S&L78 in a 
"liquid tunnel" (Schuepp and List, 1969). 

5. DISCRETE VERSION OF GAUSS' FLUX 
LAW AND SCALING RULES 

A simple method of numerically computing Sho 
using the potential V for lattice points around the 
crystal involves computation of the capacitance C. 
This is done by applying Gauss' law to compute the 

Q 

"surface charge" of the crystal, via C=— , where V 

is the surface potential assigned to the crystal at 
the onset of Jacobi's method mentioned previously. 
The total electric flux of the crystal is obtained by 
first enclosing the crystal in a "rectangular cage". 
The sides of this cage are just one lattice point away 
from the closest side of the crystal. Each side of 
this enclosure consists of its grid elements, which 
are called "area element grids". These grids have 



the same respective lattice spacing as the grid for 
the Faraday box. 

Each of the "area element" is associated with a 
set of integers /. Then the discrete approximation of 
Gauss' flux law can be written as 



E ■ dA » ~ 



dV(i,j,k) 
5ji 



■A t (11) 



which is just a sum of flux through each of the area 
elements. Note that in (11), the fact that on the 
surface of a perfect conductor, the electric field is 
always perpendicular to its surface, was utilized. The 
flux through an area element whose outward unit 
normal is in ±fi direction is approximated discretely 
by 



V((z,j,k)±fi)-V(i,j,k) 



(12) 

where v and 7 are the other two coordinates not 
equal to /1. The main reason for placing our rectan- 
gular prism enclosure at most one lattice point away 
from sides of the crystal is that such a construction 
allows dealing only with ^ rather than with the 
more complicated VV, when computing the flux. It 
is worth discussing the consequence of (12) for a spe- 
cific direction [i. For the flux through a face with its 
unit normal in ±z direction, the flux through this 
area element (12) reduces to 



F± z 



Sx 



-V(i,j,k±l) + V(i,j,k)] (13) 



The main feature of (13) is that the only unknown 
value is Sx. The flux through a face in any of the 
other directions reduce to expressions of the form 
(13). In all of these expressions, the only term with 
unknown value and dimension is 8x. For solving (8), 
it is not necessary to know the physical value of Sx 
since all the quantities, except for V, are in dimen- 
sionless form as proportionality constants a and (3. 
Thus, the numerical evaluation of potential V does 
not depend on the choice of length unit; V can be 
in whatever unit one desires. But (13) draws atten- 
tion to the actual physical value that Sx represents 
in the grid system. But note that Q ~ Sx, and so 
C ~ Sx. Thus, when computing Sh via (2), it is 



seen that Sh 



(Ml 

(5 a;) 2 



1. In other words, although 



Sx is unknown throughout our calculations, all Sx's 
cancel out at the end because similarity numbers are 



involved in our modelling. Only the relative propor- 
tions of shapes and boundary conditions are of im- 
portance. The actual physical length scales are irrel- 
evant in computation of dimensionless numbers, Sh 
and Nu, and this is well demonstrated in our model. 

The total electric flux through the rectangular 
prism enclosure is calculated by summing the flux 
(12) through each of its area elements. From the 
net flux, the "surface charge" of crystal is computed, 
and thus its capacitance. The Sherwood number 
Sh at Re=0 is computed using (2), with 5x carried 
throughout all these calculations but cancelling out 
at this last stage. 

6. RESULTS AND DISCUSSION 

The Sherwood number (Sh) for crystal shapes of 
interest in a cloud physics environment have been 
computed for zero convection (Re=0) using the 
methods outlined above. The results are listed in 
Table 1. 



Crystal 


Approximate 


Numerical 


Difference 


Shape 


Sh 


Sh 




HP 


2.74 


2.84 


+3.64% 


HC 


2.48 


2.36 


-4.83% 


BB 


4.02 


3.88 


-3.48% 


SC 


5.50 


5.70 


-3.51% 


CC 


N/A 


5.67 


N/A 



Table 1: Present numerical calculations of 
dimensionless mass transfers (Sherwood numbers) 
of hexagonal plates HP, hexagonal columns HC, 
broad branched dendrites BB, stellar crystals SC, 
and capped columns CC; comparisons with previous 
values of Sh approximated by Mc63( "Approximate 
Sh") and extrapolations by S&L78, with differences. 

The "approximate capacitances" which lead to the 
corresponding approximate Sh were given by Mc63. 
Using the analytical calculations of capacitances 
of spheres and thin circular disks of radius r to 
approximate the capacitances, Mc63, S&L78, and 
Jayaweera have given expressions for capacitance for 
the shapes listed in Table 1 (with the exception of 
CC). 



Crystal 


Approximate 


Dimensions 




capacitance 


used 


HP 


0.567r 


r=8(Sx) 


HC 




a=12.5(fe) 
b=8(6x) 


BB 


0.554r 


r=24((5.x) 


SC 


0.439r 


r=6(<5x) 



Table 2.: Approximate crystal capacitances (Mc63, 
S&L78) [cgs units] and dimensions used for the cal- 
culations in Table 1. Note that no analytical or ap- 
proximate expression of the capacitance of CC were 
given in either Mc63 or S&L78. 

In summary, the numerical methods developed and 
applied here provide convincing values of the molec- 
ular diffusion of water vapour to ice crystals of vari- 
ous shapes when compared with the values obtained 
by approximations and extrapolations (Mc63). It is 
also of interest to see that the Sherwood number for 
capped columns is not much different from the value 
for columns and that the ordering of the values for 
Re=0 is consistent with the measurements of S&L78 
at higher Reynolds numbers. 
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